Development of DGTD solver for nanophotonics applications

During the last ten years, the discontinuous Galerkin timedomain (DGTD) method has progressively emerged as a viable alternative to well established finite-difference timedomain (FDTD) and finite-element time-domain (FETD) methods for the numerical simulation of electromagnetic wave propagation problems in the time-domain. We discuss here about the development and application of such a DTGD method for solving the system of time-domain Maxwell equations coupled to material models relevant to nanophotonics. Our efforts aim at improving the accuracy, flexibility and efficiency of the method in view of the numerical treatement of realistic problems. 1. Generalities about the DGTD method The DGTD method can be considered as a finite element method where the continuity constraint at an element interface is released. While it keeps almost all the advantages of the finite element method (large spectrum of applications, complex geometries, etc.), the DGTD method has other nice properties which explain the renewed interest it gains in various domains in scientific computing: It is naturally adapted to a high order approximation of the unknown field. Moreover, one may increase the degree of the approximation in the whole mesh as easily as for spectral methods but, with a DGTD method, this can also be done very locally. In most cases, the approximation relies on a polynomial interpolation method but the method also offers the flexibility of applying local approximation strategies that best fit to the intrinsic features of the modeled physical phenomena. When the discretization in space is coupled to an explicit time integration method, the DG method leads to a block diagonal mass matrix independently of the form of the local approximation (e.g the type of polynomial interpolation). This is a striking difference with classical, continuous finite element formulations. Moreover, the mass matrix is diagonal if an orthogonal basis is chosen. It easily handles complex meshes. The grid may be a classical conforming finite element mesh, a nonconforming one or even a hybrid mesh made of various elements (tetrahedra, prisms, hexahedra, etc.). The DGTD method has been proven to work well with highly locally refined meshes. This property makes the DGTD method more suitable to the design of a hp-adaptive solution strategy (i.e. where the characteristic mesh size h and the interpolation degree p changes locally wherever it is needed). It is flexible with regards to the choice of the time stepping scheme. One may combine the discontinuous Galerkin spatial discretization with any global or local explicit time integration scheme, or even implicit, provided the resulting scheme is stable. It is naturally adapted to parallel computing. As long as an explicit time integration scheme is used, the DGTD method is easily parallelized. Moreover, the compact nature of methd is in favor of high computation to communication ratio especially when the interpolation order is increased. As in a classical finite element framework, a discontinuous Galerkin formulation relies on a weak form of the continuous problem at hand. However, due to the discontinuity of the global approximation, this variational formulation has to be defined at the element level. Then, a degree of freedom in the design of a discontinuous Galerkin scheme stems from the approximation of the boundary integral term resulting from the application of an integration by parts to the element-wise variational form. In the spirit of finite volume methods, the approximation of this boundary integral term calls for a numerical flux function which can be based on either a centered scheme or an upwind scheme, or a blend of these two schemes. In the early 2000’s, DGTD methods for time-domain electromagnetics have been first proposed by mainly three groups of researchers. One of the most significant contributions is due to Hesthaven and Warburton [1] in the form of a high order nodal DGTD method formulated on unstructured simplicial meshes. The proposed formulation is based on an upwind numerical flux, nodal basis expansions on a triangle (2D case) and a tetrahedron (3D case) and a Runge-Kutta time stepping scheme. In [2], Kakbian et al. describe a rather similar approach. More precisely, the authors develop a parallel, unstructured, high order DGTD method based on simple monomial polynomials for spatial discretization, an upwind numerical flux and a fourth-order Runge-Kutta scheme for time marching. The method has been implemented with hexahedral and tetrahedral meshes. Finally, a high order nodal DGTD method formulated on unstructured simplicial meshes has also been proposed in the same time frame by Fezoui et al. [3]. However, contrary to the DGTD methods discussed in [1] and [2], the method proposed in [3] is non-dissipative thanks to a combination of a centered numerical flux with a second-order leap-frog time stepping scheme. The extension of this DGTD method to higher-order leap-frog schemes is studied in [4]. 1.1. DGTD method for time-domain nanophotonics 1.2. Overview of related works Numerical modeling of electromagnetic wave propagation in interaction with metallic nanostructures at optical frequencies requires to solve the system of Maxwell equations coupled to appropriate models of physical dispersion in the metal. The most used are the Drude and Drude-Lorentz models. Nevertheless, although Drude and Drude Lorentz dispersion models are widely used among the numerical physicist community and proved their efficiency, they can fail to describe some materials (e.g. transition metals [5][6] and graphene [7]). Furthermore at some scales, nonlocal effects starts to play an important role [8]. When dealing numerically with Drude and Drude-Lorentz models, the FDTD method is a widely used approach for solving the resulting system of PDEs [9]. However, for nanophotonic applications, the space and time scales, in addition to the geometrical characteristics and the physical parameters of the considered nanostructures (or structured layouts of the latter), are particularly challenging for an accurate and efficient application of the FDTD method. Thus, with all their features (as described above), DG methods seem to be well suited in the context of nanophotonics. They are still not prominent compared to FDTD or FETD methods. Nevertheless, recently, unstructured mesh based methods have been developed and have demonstrated their potentialities for being considered as viable alternatives to the FDTD method [10]-[11]-[12]-[13]-[14]-[15] for nanophotonics or related applications. Besides, several studies have already been conducted regarding the development of DGTD methods for dispersive media, such as, [16]-[17]-[10]-[12]-[18] to cite a few. Noteworthy, all these studies adopt a DGTD method with upwind numerical fluxes. Furthermore one can find more studies focused on numerical analysis aspects concerning dispersive media [19]-[20]. 1.3. A non-dissipative DGTD method Towards the general aim of being able to consider concrete physical situations relevant to nanophotonics, one has to take into account in the numerical treatment, a better description of the propagation of waves in realistic media. The physical phenomenon that on has consider in the first instance here is dispersion. In the presence of an electric field the medium cannot react instantaneously and thus presents an electric polarization of the molecules or electrons that itself influences the electric displacement. In the case of a linear homogeneous isotropic media, there is a linear relation between the applied electric field and the polarization. However, above some range of frequencies (depending on the considered material), the dispersion phenomenon cannot be neglected and the relation between the polarization and the applied electric field becomes complex. In practice, this is modeled by a frequency dependent complex permittivity. Several such models for the characterization of the permittivity exist; they are established by considering the equation of motion of the electrons in the medium and making some simplifications. There are mainly two ways of handling the frequency dependent permittivity in the framework of time-domain simulations, both starting from models defined in the frequency domain. A first approach is to introduce the polarization vector as an unknown field through an auxiliary differential equation which is derived from the original model in the frequency domain by means of an inverse Fourier transform. This is called the Direct Method or Auxiliary Differential Equation (ADE) formulation. Let us note that while the new equations can be easily added to any time-domain Maxwell solver, the resulting set of differential equations is tied to the particular choice of dispersive model and will never act as a black box able to deal with other models. In the second approach, the electric field displacement is computed from the electric field through a time convolution integral and a given expression of the permittivity which formulation can be changed independently of the rest of the solver. This is called the Recursive Convolution Method (RCM). We have recently adapted the DGTD-Pp method initially introduced in [3] to deal with various dispersion models. An ADE formulation has been adopted. The resulting ADE-based DGTD method is detailed in [21] where we also study the stability and a priori convergence of the method. We first considered the case of Drude and DrudeLorentz models and, further extend the proposed ADEbased DGTD method to be able to deal with a generalized dispersion model in which we make use of a Padé approximant to fit an experimental permittivity function. The numerical treatment of such a generalized dispersion model is also presented in [21]. In this talk, we will discuss about our efforts regarding this DGTD method in order to improve its accuracy, flexibility and efficiency in view of the numerical treatement of large-scale nanophotonics applications. Numerical results will be presented for several of 3D problems ranging from academic test problems to more realistic configurations.